Dynamical density-matrix renormalization-group method 
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I present a density-matrix renormalization-group (DMRG) method for calculating dynamical 
properties and excited states in low-dimensional lattice quantum many-body systems. The method 
is based on an exact variational principle for dynamical correlation functions and the excited states 
contributing to them. This dynamical DMRG is an alternate formulation of the correction vector 
DMRG but is both simpler and more accurate. The finite-size scaling of spectral functions is dis- 
cussed and a method for analyzing the scaling of dense spectra is described. The key idea of the 
method is a size-dependent broadening of the spectrum. The dynamical DMRG and the finite-size 
scaling analysis are demonstrated on the optical conductivity of the one-dimensional Peierls-Hubbard 
model. Comparisons with analytical results show that the spectral functions of infinite systems can 
be reproduced almost exactly with these techniques. The optical conductivity of the Mott-Peierls 
insulator is investigated and it is shown that its spectrum is qualitatively different from the simple 
spectra observed in Peierls (band) insulators and one-dimensional Mott-Hubbard insulators. 

PACS numbers: 71.10.Fd, 78.20.Bh 



I. INTRODUCTION 

The density-matrix renormalization group (DMRG)Si 
is a very successful numerical methods for calculat- 
ing static properties of ground states and low-lying 
eigenstates in quantum many-body systems. For low- 
dimensional strongly correlated systems DMRG is as ac- 
curate as exact diagonalization techniques but can be 
used to study much larger systems than with those tech- 
niques (currently, up to ~ 10 3 sites). Using a finite-size 
scaling analysis it is thus possible to determine the static 
properties of a system in the thermodynamic limit with 
great accuracy. 

The calculation of dynamical properties and higher en- 
ergy excitations with DMRG has proved to be more dif- 
ficult. Several approaches have been proposed but calcu- 
lations have been carried out successfully for few prob- 
lems onb/ n The simplest of these methods is the Lanczos 
DMRG.Era In practice, this method gives accurate re- 
sults for the first few moments of a dynamical spectrum. 
Therefore, it works well for simple discrete spectra made 
of a few peaks but it usually fails for more complicated 
spectra. An alternate method for calculating dynamical 
properties is the correction vector DMRG.uu Contrary to 
the Lanczos DMRG, this method can describe complex or 
dense spectra accurately. Nevertheless, there have been 
few applicationsBQ of correction vector DMRG until now 
because this method is more difficult and requires signif- 
icantly more computer resources than the usual DMRG 
method for calculating static properties at low energy. 

In this paper I describe a simple and efficient method, 
called the dynamical DMRG (DDMRG), for calculating 
dynamical properties and excited states with DMRG. 
This approach is based on a variational principle for 
dynamical correlation functions and the related excited 
states. The variational principle is essentially an elegant 
formulation of the correction vector technique. Because 
of the variational formulation, however, the DDMRG 



method is easier to use and significantly more accurate 
than the correction vector DMRG method. 

While the spectrum of a finite system is necessarily dis- 
crete, continuous excitation bands are often found in the 
thermodynamic limit. It is possible to broaden finite-size 
spectra to simulate the continuum of an infinite-system 
spectrum. Usually, the broadening is arbitrarily large 
and no systematic quantitative analysis of finite-size ef- 
fects is performed for the spectrum. Here I show that 
dynamical properties of infinite systems can be obtained 
reliably using an appropriate finite-size scaling analysis. 
The key to the analysis is the use of a broadening which 
scales systematically with the system size. 

The DDMRG method and the finite-size scaling tech- 
nique for dynamical spectrum have already been suc- 
cessfully used to investigate the_,optical properties of 
one-dimensional Mott insulators .Ha Here I apply these 
techniques to the calculation of the optical conductivity 
in the one-dimcpaional Peierls-Hubbard model of conju- 
gated polymers£J Much effort has been devoted to under- 
standing the optical properties of these materials. In par- 
ticular, the optical conductivity of the Peierls-Hubbard 
model has been studied extensively and analytical results 
have been obtained for various special cases, such as the 
Mott-Hubbard insulator and the Peierls (band) insula- 
tor limits. Leaving out these special limits the Peierls- 
Hubbard model describes a Mott-Pcicrls insulator and 
its optical properties are still poorly understood. In this 
paper I show that DDMRG can accurately reproduce the 
known analytical results for the optical spectrum in the 
thermodynamic limit. Then I investigate the optical con- 
ductivity of a Mott-Peierls insulator using DDMRG and 
show that it displays specific features. 

The paper is organized as follows: The variational prin- 
ciple for dynamical correlation functions and related ex- 
cited states is presented in the next section. I describe 
the dynamical DMRG method in Sec. E. The finite- 
size scaling analysis is presented in Sec. [Vl I report 
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and discuss the results for the optical conductivity of the 
Peierls-Hubbard model in Sec. |V|. Finally, I conclude in 
the last section. 



II. VARIATIONAL PRINCIPLE 

The dynamic response of a quantum system to a time- 
dependent perturbation is often given by dynamical cor- 
relation functions such as 

G a (oj + ir,) = —{fh\Al - 1 . -A\%) , (1) 

7T Eq + U! + tT] — H 

where H is the time- independent Hamiltonian of the sys- 
tem, E and |^q) are its ground-state energy and wave- 
function, A is the quantum operator corresponding to 
the physical quantity which is analyzed, and A^ is the 
Hermitian conjugate of A. A small real number r\ > 
is used in the calculation to shift the poles of the corre- 
lation function into the complex plane. (I set U = 1 in 
sections II to IV.) 

In general, we are interested in calculating the imagi- 
nary part of the correlation function 



I A (w + irj) = Im Ga + irj) 
1 



(2) 



n 



;A\^o) 



in the limit r\ — > 



I A (ui) = lim I A (u> + > 

7)^0 



(3) 



It should be noted that the spectrum I a (u> + irj) for any 
finite r\ > can be calculated from the spectrum I A {u>) 
by convolution with a Lorentzian distribution 



I A {uj + iri) = C v [I A {w)) > , 



(4) 



where I use the notation C v [f(w)] to represent the con- 
volution of a spectral function f(u>) with a Lorentzian 
distribution of width r] 



c v [m] = 



dw'f(u>')- 



n 



(5) 



7T (lu — uj') 2 + r/ 2 
The moments of the spectrum Ia(w) fulfill sum rules such 



— oo 



duI A (u) = {ik\A*A\1>o) , 

duI A (u)w = {tfo\A*[H,A]\fo) , 

+oo 

dujI A {u)uj 2 = (M[A\H][H,A]\i, ) , 



(6) 



A dynamical correlation function (nj) can be calculated 
using the correction vector method. The correction vec- 
tor associated with Ga (w + in) is defined by 



\ip A (u; +irj)) 



1 



E + u) + irj — H 



\A) 



(7) 



where \ A) — A\ipo). If the correction vector is known, the 
dynamical correlation function can be calculated directly 



G a (uj + in) = --(A\i/j a (uj + itj)) 

7T 



(8) 



To calculate a correction vector one first solves an inho- 
mogeneous linear equation 



[(E + uj - H) 2 + if] \$) = -n\A) 



(9) 



which always has a unique solution \ip) = |Ya(u; + in)) 
for rj 7^ 0. The correction vector is then given by 

\ip A (oj + in)) = \x a (oj + in)) + i\Y A (uj + in)) (io) 

with 

\x A (u + i v )) = H ~ E ° ~ " \y a (u + in)) ■ (ii) 
n 

One should note that the states \X A (ui+in)) and 
in)) are complex if the state \A) is not real, but they 
always determine the real part and imaginary part of the 
dynamical correlation function Ga(w + in), respectively, 



BbG a (u + %v) = --(A\X A (oj + in)) , (12a) 
ImG A {uj + in) = (A\Y a (uj + in)) . (12b) 

The derivatives of the real and imaginary parts can also 
be calculated from these states 

— ReG^ + w?) = -[{X A {u + in)\X A {u + in)) 
aoj ir 

-(Y A (oJ + in)\Y A (u + in))} (13) 

— ImG^C^ + in) = —(X a (uj + in)\Y A (u; + in)) ■ 
auj ir 

A well-established approach for solving an inhomoge- 
neous linear equation (|9|) is to formulate it as a minimiza- 
tion problem. One considers the functional 

W A<V M) = (i>\(E +co-H) 2 +n 2 W 

+n{A\ip) + n{ip\A) . (14) 

For any n ^ and a fixed frequency uj this functional 
has a well-defined and non-degenerate minimum for the 
quantum state which is solution of Eq. (^) 

(15) 



|V>min) = \Ya(u + in)) 



It is easy to show that the value of the minimum is 
related to the imaginary part of the dynamical correlation 
function 



where [A, B] = AB - BA. 



W a , v (uj, ipmin) = -nr)I A (w + irj). 



(16) 
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Therefore, the calculation of spectral functions can be 
formulated as a minimization problem. To determine 
7,4 (a; +177) at any frequency u and for any 77 > 0, one min- 
imizes the corresponding functional W A ,ri{w, ip)- Once 
this minimization has been carried out, the real part of 
the correlation func tion Ga(w + irj) can be calculated 
using Eqs. ( |TT| ) and (12a) if necessary. This is the varia- 
tional principle for dynamical correlation functions. 

It is clear that if we can calculate \Y A (u) + ir])) exactly, 
this variational formulation is completely equivalent to 
the correction vector method. However, if we can only 
calculate an approximate solution with an error of the or- 
der e < 1, = \Y A {u + irj)) + e\<j>) with {<j)\<j)) = 1, the 
variational formulation is more accurate. In the correc- 
tion vector method t he er ror in the spectrum I a (oj + irj) 
calculated with Eq. (12b) is also of the order of e. In 
the variational approach it is easy to show that the error 
in the value of the minimum Wa^m, V'min), and thus in 
Ia(w + iff), is of the order of e 2 . With both methods the 
error in the real part of Ga{^> + iff) is of the order of e. 

One can write the function I A {w) in the spectral form 
(or Lehmann representation) 

I A (u>) = ^|<^ n |^o>| 2 <$(u; + E a - En) , (17) 

n 

where |?/>o) is the ground state, \ip n ),n > 1, denotes the 
other eigenstates of H, and E 0} E n are their respective 
eigenenergies. Obviously, only the eigenstates with a 
finite matrix element (^ n |^4|^o) 7^ contribute to the 
spectrum and here we are only interested in those excited 
states. In the correction vector method the excitation en- 
ergies E n — Eq and the spectral weights KVvJ^IV'o)! 2 can 
be obtained from the poles of Ga(u + iff)- The corre- 
sponding wavefunctions |^ n ) can be calculated by taking 
the r\ — ► limit of the correction vectors 



\ip n ) cx lim \Y A {E n 

rj—>0 



Ea + if})) • 



(18) 



The excited states contributing to G A (^> + iff) corre- 
spond to the local maxima of the spectrum Ia{u + iff) 
for small enough r\ > 0. Therefore, they can also be ob- 
tained by minimization of the functional Wa.ti(^, VO with 
respect to both to and ip. The local minima of Wa,t)(u, tji) 
are given by the conditions 



Eq 



(VwlVw) 



(19) 



In the limit 77 — > 0, w m in + Eq tends to the energy E n of 
an eigenstate with finite spectral weight, |f/ , min) is equal 
to the corresponding eigenstate \ip n ) up to a normaliza- 
tion constant, and — WA^Wmin, V'min) tends to the spec- 
tral weight I {ijj n \ A\%j]Q) | 2 . This is the variational principle 
for excited states contributing to a dynamical correlation 
function Ga {oj + iff) . 

Again this variational formulation is completely equiv- 
alent to the correction vector method if |Ya(w + if})) can 



be calculated exactly. In an approximate calculation, 
however, errors in the eigenenergies and spectral weights 
are of the order of e with the correction vector method, 
while they are of the order of e 2 with the variational for- 
mulation. 



III. DYNAMICAL DMRG 

DMRG is a numerical method for calculating the prop- 
erties of lattice quantum many-body systems. It is de- 
scribed in detail in several publications (for instance, see 
Refs.fi] and |). DMRG can be considered as a variational 
approach. The system energy 



(20) 



is minimized in a variational subspace (the DMRG ba- 
sis) of the system Hilbert space to find the ground-state 
wavefunction \ipo) and energy E Q — E(ip Q ). If the ground- 
state wavefunction is calculated with an error of the order 
of e < 1 (i.e., \ip) = IV'o) + e\4>) with (</>]</>) = 1), the en- 
ergy obtained is an upper bound to the exact result and 
the error in the energy is of the order of e 2 (as in all 
variational approaches). In principle, the DMRG energy 
error is proportional to the weight of the density-matrix 
eigenstates discarded in the renormalization procedure. 
This discarded weight can be reduced by increasing the 
number to of density-matrix eigenstates kept in the cal- 
culation, which corresponds to an increase of the varia- 
tional subspace dimension. Therefore, the energy error 
systematically decreases with increasing to in a DMRG 
calculation. 

The DMRG procedure used to minimize the energy 
functional ( pp[ ) can also be used to minimize the func- 
tional W AjTI (uj,if>) and thus to calculate the dynamical 
correlation function G a (uj + iff). I call this approach 
the dynamical DMRG method. The minimization of the 
functional Wa,ti(^^) is easily integrated into the usual 
DMRG algorithm. At every step of a DMRG sweep 
through the system lattice, a superblock representing the 
system is built and the following calculations are per- 
formed in the the superblock subspace: 

1. The energy functional E(ip) is minimized using 
a standard iterative algorithm for the eigenvalue 
problem. This yields the ground state vector \ip ) 
and its energy E in the superblock subspace. 

2. The state \A) is calculated. 

3. The functional tp) is minimized using an 
iterative minimization algorithm. This gives the 
first part of the correction vector |Ya(w + irj)) and 
the imaginary part Ia{w + iff) of the dynamical 
correlation function. 

4. The second part \X A (^ + if])) of the correction vec- 



tor is calculated using Eq. (TTTl) 
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5. The real part and the derivatives of the dynam- 
ical correlation function can be calculated from 
Eqs. (12a) and (|l3|), respectively. 



6. The four states \ip ), \A), \Y A (oj+ir))), and \X a (uj + 
irj)) are included as target in the density-matrix 
renormalization to build a new superblock at the 
next step. 

The robust finite-system DMRG algorithm must be used 
to perform several sweeps through a lattice of fixed size. 
Sweeps arc repeated until the procedure has converged to 
the minimum of both functionals E(ip) and Wa. v {uj, ip). 

To obtain the dynamical correlation function G a {oj + 
irj) over a range of frequencies, one has to repeat this cal- 
culation for several frequencies ui. If the DDMRG calcu- 
lations are performed independently, the computational 
effort is roughly proportional to the number of frequen- 
cies. It is also possible to carry out a DDMRG calculation 
for several frequencies simultaneously, including several 
states \X A (ui + ir/)) and |Ya(^ + irj)) with different fre- 
quencies u) as target. The optimal number of different 
frequencies to be included in a single calculation depends 
strongly on the problem studied and the computer used. 
As calculations for different frequencies are essentially in- 
dependent, it would be easy and very efficient to perform 
these calculations on a parallel computer. 

If one performs a DDMRG calculation for two close 
frequencies u>% and u>2 simultaneously, it is possible to cal- 
culate the dynamical correlation function for additional 
frequencies lo between wi and 102 without including the 
corresponding states \XA(u>+irf)} and \Y A (u)+irj)} as tar- 
get in the density-matrix renormalization. This approach 
can significantly reduce the computer time necessary to 
determine the spectrum over a frequency range but the 
results obtained for w ^ u>x,u>2 are less accurate and not 
always reliable, as the DMRG basis is optimized for the 
frequencies io\ and u>2 only. (A similar technique is the 
calculation of spectra with the Lanczos algorithm in the 
DMRG basis optimized for a pair of correction vectors, 
see Ref. ||.) Alternatively, between the frequencies for 
which Ga{^ + iv) is determined directly with DDMRG, 
we can calculate the dynamical correlation function by 
interpolation using the DDMRG data for the function 
and its derivative. 

If a complete spectrum Ia(<-o + if)) has been obtained, 
it is possible to calculate the moments of the spectral 
distribution [the left-hand-side of Eq. (6)]. The first few 
moments [the right-hand-side of Eq. (6)] can be calcu=. 
lated accurately using the Lanczos DMRG methodE'cl 
This provides an independent check of DDMRG results. 
[Note that only the first sum rule (||) is satisfied exactly 
for 77 > 0.] 

To calculate individual excited states contributing to 
the spectrum in a given frequency range (u>i, L02), one in- 
cludes a minimization of Wa^I^, t/>) with respect to 10 
(u>i < u> < L02) in the third step of the DDMRG algo- 
rithm described above. In this case, |YA(w m m + and 
|^(w m in + ill)) are included as target in the sixth step. 



The parameter 77 must be much smaller than the dis- 
tance E n+ i — E n between two successive eigenstates con- 
tributing to the dynamical correlation function or must 
decrease during the calculation until the desired accu- 
racy is obtained. To make the procedure robust it is 
necessary to simultaneously target a second correction 
vector \ip A (u + i?])) with a fixed frequency and a param- 
eter 77 of the order of the frequency range. Typically, I 
use to = (ll>i + u>2)/2 and 77 = (u>2 — Wi)/4. 

Because of the variational principle one naively ex- 
pects that the DDMRG results for I a + irj) must con- 
verge monotonically from below to the exact result as the 
number m of density-matrix eigenstates is increased. In 
practice, the convergence is not always regular because 
of two approximations made to calculate the functional 
WAtf{oj,ilS) in a DMRG basis. First, the ground-state 
energy Eq and the state \A) used in the definition (|l4| ) 
of Wk,tj(w, VO are n °t known exactly but calculated with 
DMRG. If the number m of density matrix eigenstates 
is increased, Eq and \A) are modified (they become pro- 
gressively more accurate) and the functional Wa, v (u,iP) 
is changed, which can affect its minimum arbitrarily. We 
also note that errors of the order of e in Eq or \A) re- 
sult in errors of the same order in Ia(w + iv)- Therefore, 
to observe a regular convergence with increasing m and 
to obtain accurate results for I a (w + 777) , it is necessary 
in the first place to determine the ground state and the 
state I A) with great precision (and thus to include the 
state I A) as target). 

To calculate the functional Wa^(w, ijj) in the third step 
of the DDMRG algorithm, one needs an effective repre- 
sentation of the operator (H — Eq — lo) 2 in the superblock 
subspace 

[(H - E Q ~ cjfU = 0\H -E -lv) 2 O , (21) 

where the operator O represents the projection onto the 
superblock subspace. For a typical many-body Hamil- 
tonian H such a calculation is excessively complicated 
and computationally intensive. Therefore, I calculate an 
effective representation of H only, H c g — O^HO, and 
assume that 



[{H -Eq- tof 



eff 



(-ffoff — e — 



(22) 



to calculate Wa, v (uj, ip) in the superblock subspace. This 
second approximation can cause a violation of the varia- 
tional bound W A v {u), ip) > — 7^7/^4 (u + 777). Fortunately, 
the substitution (E2p has no significant effect on the mim- 
mumof Wa^^,^) if the state (H— Eq— ui)\YA(oj+ir))) oc 
\XA(oj+irj)} is accurately represented in the DMRG basis 
[i.e., if 01X^4(0; + 777)) w + 777)) for all superblock 

sub-spaces]. Therefore, to use the substitution ( p2] ) with- 
out loss of accuracy it is necessary and sufficient to in- 
clude the state + 777)) as a target in a DDMRG 
calculation, even if the real part of the dynamical corre- 
lation function is not calculated. 

In practice, for sufficiently large m I have found that 
the absolute values of errors in a spectrum Ia{^> + iff) 
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decrease systematically with increasing m. Therefore, it 
is possible to estimate the accuracy of a DDMRG cal- 
culation from the results obtained for different values of 
m as one can do for static properties calculated with 
DMRG. Moreover, DDMRG results for I a (uj + it]) tend 
to be smaller than the exact result for almost all frequen- 
cies although they can exceed it occasionally. 

Obviously, DDMRG is very similar to the correction 
vector DMRG.beI The same DMRG basis is built because 
the same target states are used in both methods. As nu- 
merical errors are usually dominated by the DMRG basis 
truncation, the correction vector parts + irf)) and 

|Ya(w + irj)} are calculated with the same precision e in 
both methods for a given number m of density-matrix 
eigenstates kept per block. Nevertheless, the variational 
formulation has two significant advantages. First, the 
errors in the spectrum I a (to + it]), the excitation ener- 
gies E n — Eo, and the spectral weights are of the order 
of e 2 instead of e in the correction vector method, as 
explained in Sec. II. If one uses the Lanczos algorithm 
instead of Eq. (||) in the correction vector DMRG, errors 
become even larger. Therefore, DDMRG results are more 
accurate than those obtained with the correction vector 
DMRG for a given number m of density-matrix eigen- 
states. Second, a DDMRG calculation is essentially an 
application of the standard DMRG algorithm to the min- 
imization of a different functional. In particular, the nu- 
merical accuracy and computational effort are controlled 
by the sole parameter m in an optimized DDMRG cal- 
culation as in a ground-state DMRG calculation_|-|Thc 
correction vector DMRGo'El and Lanczos DMRGeIq are 
significantly more complicated than the standard DMRG 
method. In particular, the numerical accuracy and com- 
putational effort depend significantly and sometimes un- 
predictably on the specific states (Lanczos vectors and 
correction vectors) included as target in the density- 
matrix renormalization. Therefore, it is easier to imple- 
ment and use DDMRG than the correction vector DMRG 
or the Lanczos DMRG. 



IV. FINITE-SIZE SCALING 

DDMRG allows us to calculate spectral functions of 
a correlated electron (or spin) system on a finite lattice 
with a broadening given by the parameter 77 > 0. They 
have the generic form 

^M = ^..W kW ^ + ^ (23) 

where Lo n (N) denotes the excitation energy and A n (N) > 
the spectral weight of the system eigenstates, and TV is 
the number of lattice sites. Such spectra are discrete for 
77 — ^ because there is only a finite number of eigenstates 
in a finite system. In the thermodynamic limit a spectral 
function 

I (to) = lim lim In v (u)) (24) 

j;^0 JV— >oc 



can contain discrete and continuous parts. (It should be 
noted that the order of limits in the above formula is 
important.) To determine the properties of a dynamical 
spectrum I(uj) in the thermodynamic limit one has to 
analyze the scaling of the corresponding spectra //v,r;(w) 
as a function of system size. Here I present a finite-size 
scaling technique for spectral functions calculated with a 
numerical method such as DDMRG. 

Computing both limits in Eq. ( pl| ) from numerical re- 
sults for Ipf tV (uj) requires a lot of accurate data for dif- 
ferent values of ij and iV and can be the source of large 
extrapolation errors. A much better approach is to use a 
broadening r](N) > which decreases with increasing TV 
and vanishes in the thermodynamic limit. The dynamical 
spectrum is then given by 

I(u) = lim I nmn) (uj) (25) 

N — too 

= ^ \ E A " W ( LJn ( N ) 11 %Kv 2 (N) ■ 

From the existence of both limits in Eq. ( p4[ ) it can 
be demonstrated that there exists a minimal broaden- 
ing 770 (TV) > 0, which decreases as a function of TV and 
converges to zero for TV — > 00, such that the above equa- 
tion is valid for all functions rj(N) with rj(N) > 770 (TV) 
and limjv-,00 v(N) = 0- The function 770 (TV) depends on 
the frequency u> considered. For a finite lattice with TV 
sites, let M Wj£ (TV) be the number of excited states con- 
tributing to the spectral function in a small interval of 
width e around the frequency to (i.e., \uj n (N) — uj\ < e/2). 
If M Wj£ (TV) remains finite for any e > as TV — > 00, the 
spectral function I(lo) is discrete at u> and 770 (TV) = 0. 
Equivalently, one can take the 77 — i> limit first in 
Eq. Q. If M u , e (iV) diverges for all e > as TV -> 00, 
the spectrum is dense at lu and a minimal broaden- 
ing 770 (TV) > is required for Eq. (p5|) to be valid. 
For instance, 770 (TV) must be larger than the distance 
Slu = u n+ i(N) — uj n (N) between two consecutive excited 
states in the spectrum. Note that while a continuous 
spectrum is obviously dense, a dense spectrum can be 
continuous or discrete. For instance, an infinite number 
of excited states with A n (N) > can converge to the 
same energy as N — > 00. This seems to happen for the 
optical conductivity associated with an exciton in a open 
chain .□ 

The function f?o(TV) depends naturally on the specific 
problem studied [i.e., the scaling of the energies U) n (N) 
and spectral weights A„(TV)]. For the optical conductiv- 
ity of one-dimensional correlated electron systems such 
as the Peierls-Hubbard model, I have found numerically 
that a sufficient condition for all frequencies w in a dense 
part of the optical spectrum is 



where the constant c is comparable to the width of the 
dynamical spectrum I(u), which is finite in such lattice 
models. Usually, one wants to keep the broadening r\ as 
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small as possible because it reduces the resolution of the 
spectrum. Therefore, I use 

V(N) = £ (27) 

in Eq. ( p5| ) to analyze the finite-size scaling of spectral 
functions calculated with DDMRG and to extrapolate 
the finite-size results to the thermodynamic limit. The 
condition (^6|) has a very simple physical interpretation. 
The spectral function In^uj) represents the dynamical 
response of the system over a time period ~ I/77 after 
one has started to apply an external force. Typically 
in a lattice model the spectral width is proportional to 
the velocity of the excitations involved in the system re- 
sponse. Thus the condition ( p6| ) means that excitations 
are too slow to travel the full length ~ N of the system 
in the time interval ~ \jr\ and do not "sense" that the 
system is finite. 

An additional benefit of a broadening satisfying the 
condition (pif ) is that the finite-system spectrum In jV (u>) 
becomes indistinguishable from the infinite-system spec- 
trum with the same broadening 77 for relatively small N, 

I N , v (u) » C„[J(w)]. (28) 

Therefore, if one knows an exact or conjectured spec- 
tral function I{ui) for an infinite system, its convolution 
with a Lorentzian of width 77 can be compared directly 
with the numerical results for the finite-system spectrum 

In. 

r){uj). This approach has been applied successfully to 
the optical conductivity of one-dimensional Mott insula- 
tors in Refs. |^ and |^, and additional examples are pre- 
sented in the next section. 

In practice, the extrapolation scheme ( p5j ) works well 
at fixed frequency for the continuous parts and the non- 
dense discrete parts of a spectrum I(ui) only. To detect 
singularities in I(oj) and determine their properties, it 
is generally easier to analyze the scaling of maxima in 
In. 

v (ui) or in its derivative as N — > 00 and 77 — ^ 0. To 
perform this scaling analysis one can use a size-dependent 
broadening 7]{N) such that Eq. (53) is valid and In, v (uj) 
is a good approximation of C v [I(uj)] around the maxi- 
mum. Then the scaling of a maximum in //v )r) (w) for 
rj(N) — > gives the scaling of the corresponding maxi- 
mum in Cjj [I(u>)] for 77 — > 0. Here I discuss some examples 
of this technique which arc useful for the analysis of the 
Peierls-Hubbard model optical conductivity presented in 
the next section. [One should also note that to detect the 
presence of a gap between two bands in an infinite-system 
spectrum I(u>), it is often faster and more reliable to in- 
vestigate the scaling of minima in In jTI (u)) for 77(A) — > 
than to perform extrapolations at fixed frequencies.] 

First, we consider an infinite-system spectrum with a 
peak in a continuous band, 

I{i0) = I Q 6(0J -L0q)+ /cent (w) (29) 

for \lo— ujq\ < A, where Iq > and I CO nt(u) is a continuous 
function. It is easy to show that for 77 <C A the maximum 



of C v [I(ui)) diverges as Iq/{ttt]) and that the position of 
the maximum converges to luq for r\ — > 0. The maximum 
of the corresponding finite-system spectra has 
the same scaling properties for r/(N) — > 0. Therefore, 
it is possible to detect such a <5-peak and determine its 
weight Iq in an infinite-system spectrum, even if Iq is 
only a small fraction of the total spectral weight of I CO nt ■ 
Second, we consider an infinite-system spectrum with 
a power-law divergence at the band edge 

J(w) = I 6(u - wo)|w - u \- a (30) 

for \u> — luq\ < A, where 9{x) = for x < and 8(x) = 1 
for x > 0, Iq > 0, and < a < 1. One can show 
that for r\ <C A the maximum of C v [I(uj)] diverges as 
?y _Q and that the position of the maximum converges to 
u>o from above for 77 — > 0. Again, the maximum of the 
corresponding finite-system spectra 7jv,jj(w) has the same 
scaling properties for r)(N) — * 0. Thus it is possible to 
detect such a singularity and determine the exponent a 
from the finite-system numerical data. 

Third, we consider a continuous infinite-system spec- 
trum with a singularity in its derivative 

I(iu) = I 6{u;-u;o)\u;-LUo\ a (31) 

for \lu — ujq\ < A, where the function 0(x) and the con- 
stants Iq and a are as in the previous example. Here the 
derivative of C v [I(lo)] has a maximum which diverges 
as rj 01 ^ 1 for rj -C A while its position converges to ujq 
from above as 77 — > 0. The maximum in the derivative 
of /jv^(w) has the same scaling properties for 77^) — » 0. 
Therefore, it is possible to determine the exponent a from 
the finite-system numerical data in this case too. 

Finally, we consider a special function representing a 
continuous spectrum above a gap and a truncated diver- 
gence close to the band edge 

j ( \ T al 2y/uj Q \uj - UJ \ 

I(w) = Iq0(uj - u ) — (32) 

7W0 + — too] 

for \ui — wo I < A, where 7 is a constant such that 
< 7 < A/ujq. The function 9{x) and the other con- 
stants are as in the previous examples. For 7 > this 
spectrum vanishes as a square- root at the band edge luq, 
goes trough a maximum Iq j^Ff at u> — (1 +7)0*0 then de- 
creases monotonically. For 7 <C 1 the maximum is very 
sharp and close to the band edge, and I(lu) appears to di- 
verge as 1 / -y/jaJ^tJoI at higher energy uj > (1+7)^0 • The 
continuous vanishing of I(u>) at the band edge is appar- 
ent only in a small frequency range u>q < u> < (1 + 7)0^0 ■ 
As 7 — > this maximum becomes a square-root singular- 
ity at ujq. A qualitatively similar behavior is often found 
in the optical conductivity of one-dimensional insulators 
(see the discussion in the next section). Obviously, the 
maximum of C v [I(lo)\ tends to To/i/7 an( i its position 
converges to oj = (1 + 7)^0 for 77 — * 0. For 7 < 1, 
however, the convergence of the maximum becomes ap- 
parent only for 77 <C 7^0- For larger 77 the maximum 
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appears to diverge as l/y 7 ??. Similarly, the maximum of 
the derivative diverges as 1/y/rj for 77 —* as discussed 
the previous example. In the present case, however, this 
scaling is not observed as soon as 77 <C A but only if 
■q <C 7cjo- The finite-system spectrum I^ jn (uj) and its 
derivative have the same scaling properties for r)(N) — > 0. 
Therefore, with the scaling analysis of finite-system spec- 
tra it is possible to distinguish an infinite-system spec- 
trum with a truncated divergence above the band edge 
from a spectrum with a real divergence at the band edge, 
provided that one can do calculations with a resolution 
r](N) < 7^0 • 

In summary, the dynamical spectrum of an infinite sys- 
tem can be determined accurately and efficiently from 
numerical data for finite-system spectra using a size- 
dependent broadening rj(N). The broadening r)(N) must 
be larger than a minimal broadening r]o{N), which de- 
pends on the system investigated and can vary with the 
frequency. Often this broadening conceals the finite-size 
effects and one can directly compare finite-system spec- 
tra to analytical results for infinite systems using a con- 
volution with a Lorentzian distribution, see Eq. (p8[). If 
this comparison is not possible or not sufficient, specific 
points of the spectrum can be extrapolated to the ther- 
modynamic limit using Eq. (25). Finally, the scaling of 
maxima in finite-system spectra or their derivatives [for 
i](N) — > 0] allows us to find and analyze singularities in 
the infinite-system spectrum. For one-dimensional corre- 
lated electron systems a sufficient condition for the min- 
imal broadening is given by Eq. (pfjf) and one can use a 
size-dependent broadening (|27 



V. OPTICAL CONDUCTIVITY OF THE 
PEIERLS-HUBBARD MODEL 

In this section I apply the DDMRG method and the 
finite-size scaling analysis to the optical concbictivity 
of the one-dimensional Peierls-Hubbard model.Eij This 
model is defined by the Hamiltonian 



H = T 



1=1 



(33) 



with 



T 



( 34 ) 

It describes electrons with spin a =|, | which can hop be- 
tween neighboring sites in a lattice with an even number 
N of sites. In Eq. ( |34| ) the index I runs from 1 to N — 1 
for an open chain and from 1 to N if periodic bound- 
ary conditions are used. Here c~l a , ci^ a are creation and 
annihilation operators for electrons with spin a at site I 
and n^ a — cf a ci, a are the corresponding density opera- 
tors. The hopping integral t > gives rise to a single- 
electron band of width At. The dimerization parameter 



< |A| < 2t determines the strength of the periodic lat- 
tice potential generated by the Peierls instability. (For 
a finite open chain I only use A > to avoid spurious 
excitations at the chain ends.) The Coulomb repulsion 
is mimicked by a local Hubbard interaction U > 0. The 
number of electrons equals the number of lattice sites. 

The ground state, single-particle charge gap and spin 
gap of this system can be calculated with great accuracy 
on lattices with up to N ~ 10 3 sites using DMRG. The 
single-particle charge gap is given by 

E C {N) = E (N + 1) + E (N - 1) - 2E (N) , (35) 

where Eq(M) denotes the ground-state energy for M 
electrons in a -/V-site system. For an even number N 
of sites and electems the Peierls-Hubbard model ground 
state is a singletliil and the spin gap is given by 



E S (N) = E {S Z = ft) - E {S Z = 0) 



(36) 



where S z is the the z-component of the total spin and 
E (S Z ) is the ground-state energy for a fixed value of S z . 

Spectroscopy with electromagnetic radiation is a com- 
mon experimental probe of solid-state materials. The 
linear optical absorption is proportional to the real part 
(Ti(u)) of the optical conductivity. For to > 0, <Ji(u)) is 
related to the imaginary part of the current-current cor- 
relation function Gj(huj + irj) by 



<7i(u>) 



lim Im Gj(hu> + irj) 



Nau 17-1-0 

^-Y,\^\J\^n)\ 2 6(hLJ + E -E n ) 



(37) 



Here \i/jq) is the ground state of the Hamiltonian H, 
\ip n )(n > 1) are the other eigenstates of H, and E , E n 
are their respective eigenenergies. In this model the cur- 
rent operator is 



J 



lae 



E *-(-!)■ 



A 



c 



+ 

1+1,(7 



CL, 



(38) 

where a is the lattice constant, — e is the charge of an 
electron, and the index I takes the same values depending 
on the boundary conditions as in Eq. (|34|). Note that this 
is the natural definition of the current operator for both 
types of boundary conditions. The Parzen filter used 
for open boundary conditions in other workstiQ is not 
necessary and thus not used in this work. 

In an open chain the optical absorption is also related 
to the dynamical polarizability a(u>), which is given by 
the imaginary part of the dipole-dipole correlation func- 
tion Gd(Tiu + irj), 



a(u) = lim Im Gd(^ + iff) 

Na ?j->o 



(39) 



?-J2\(iPo\D\iP n )\ 2 5(hu> + E -E n ) 



Na 
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For a constant lattice spacing a the dipole operator is 

N 



D = -eaJ^Hni - 1) 



1=1 



with ri; = n/j + ni t i. Using the relation J 
one easily proves that 



(40) 



-mm 



o"i (u>) — uia(ui) 



(41) 



and 



crUu) = lim Im{(Ku + in) G D (huj + in)} . (42) 

Nan 7j->o 

The optical conductivity can be calculated with the 
DDMRG method described in this paper. For an open 
chain Eqs. (37),(|39|), and ( |4^ ) provide us with three dif- 
ferent approaches. First, one can calculate the imagi- 
nary part of the current-current correlation function with 
DDMRG and use Eq. (|37]) to obtain the convolution of 
the reduced optical conductivity 



C v [wax(u)} 



7T 

Na 

-y— 

Na ^ (ha; 



Im Gj(huj + irj) 

r?|(V>o|J|^)| 2 
+ E - E n y + rj 2 



(43) 



This is also the only approach possible with periodic 
boundary conditions. Second, one can calculate the 
imaginary part of the dipole-dipole correlation function 
with DDMRG and use Eq. (|39| ) to obtain the convolution 
of the dynamical polarizability 



C,[aM] 



Im GdOt-u 

Na 



in) 



(44) 



Na ^ (Huj + E - E n ) 2 



T) 



The optical conductivity is then given by the rela- 
tion (|l|). Finally, one can calculate the complete 
dipole-dipole correlation function with DDMRG and use 
Eq. (E2) to obtain the convolution of <ti(lu) directly 



C n [ai (oj)] 



Nafi 
1 

Nafi 



Im{(huj + in) G D {huj + ii])} (45) 

n\(^\D\^ n )\ 2 {E n -E Q ) 



(Tioj + E Q - E n ) 2 



vr 



C n [ai(uj)] can also be formulated in terms of the current 
matrix elements | (ipo\J\tp n ) \ 2 [see Eq. (4) of Ref. ||. Al- 
though the real part of a dynamical correlation function 
is used in Eq. ( f45| ) to calculate the optical conductiv- 
ity, its relative contribution to C n [<Ji{ijj)\ is of the order 
(ij/t) 2 . Therefore, the numerical precision is not signifi- 
cantly reduced by the lower accuracy of DDMRG for the 
real part of dynamical correlation functions. 

Clearly, all three approaches give the same spectrum 
a\(uo) for 77 — > 0. In DDMRG calculations with 77 > 0, 



however, they are not equivalent. First, I have found 
that it is easier to calculate the dipole-dipole correlation 
function Go{nLJ + irf) than the current-current correla- 
tion function Gj{Tioj + irj), except for very strong cou- 
pling U 3> t. Second, the finite-size scaling [using a size- 
dependent broadening (^)] is different for the three op- 
tical spectra C v [uj(Ji(uj)]/u), wC,[a(w)[, and C v [ai(uj)], 
especially for very small and very large frequencies u. 
Usually, C v [ai(uj)] is the best approximation to <7\(lo) 
but at low energy {Tiuj < t) it can be more convenient 
to use wC,[a(w)] while at high energy (Huj ^> t) I prefer 
C v [uai(oj)}. 

To calculate the optical spectrum of the Peierls- 
Hubbard model I have used the third approach, Eq. (|4"s|), 
in most cases. Thus C n [cri(uj)] is shown in the figures of 
this paper unless I state explicitly otherwise. Only op- 
tical spectra calculated with open boundary conditions 
are presented. As the size-dependent broadening ( |27| ) 
conceals most of the finite-size effects, spectra calcu- 
lated with periodic boundary conditions would be al- 
most identical. In all figures showing optical spectra I 
set a — e — h = t = \. Thus <j\{uj) is shown in units of 
e 2 a/h, Loai{uj) in units of e 2 at/h 2 , and the frequency to 
in units ofU/t. 

The sum rules (^) take a simple form for the opti- 
cal conductivity in the Peierls-Hubbard model with open 
boundary conditions 



dujuj(ji{u)) 



du)<Ji(u>) = 


dco^l = 
^ 



=^ <*,m*,>, 



(46a) 
(46b) 
(46c) 



To prove the second sum rule (46b) one uses the rela- 
tion [D,J] = —ia 2 e 2 T/h. With periodic boundary con- 
ditions, only the first two sum rules remain valid. In the 
second sum rule, however, one must take into account 
the coherent part of the conductivity at to — and the 
proof is more complicated than for an open chain. 

The right-hand-side of Eq. (^6|) can be calculated accu- 
rately with the ground-state or Lanczos DMRG method. 
Note that for 77 > an optical spectrum calculated from 
Eq. @, (Es[K or @) exactly fulfills the sum rule ( |46a] ) , 
( p6h| ) , or p6c] ), respectively. For the DDMRG results 
presented in this paper, the sum rules are fulfilled within 
a few percents. 

The optical gap E op t(N) is defined as the excitation 
energy (E n — Eq) of the lowest eigenstate contributing 
to the optical conductivity (i.e., (V'nl^lV'o) 7^ 0) in a N- 
site system. E opt (N) can be calculated with the DDMRG 



method for individual excited states described in Sec. Ill 
As the Peierls-Hubbard Hamiltonian (^) has a particle- 
hole symmetry the optical gap can alsoJpe determined 
using the symmetrized DMRG methodic As expected, 
both approaches give the same results for E opt (N) within 
numerical errors. In the thermodynamic limit (N — > 00) I 
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have found that the optical gap E opt is equal to the single- 
particle charge gap E c (||) for all U > and 2i > A > 0. 
[In the dimer limit (A = 2t) the Hamiltonian ([33]) de- 
scribes independent dimers: the optical weight is con- 
centrated in a single peak corresponding to Frenkel ex- 
citons localized on a dimer and the other (delocalized) 
excitations above and below this peak carry not optical 
weight, thus E c < -E^pt-] In a-finite system or with addi- 
tional electronic interactions^ the single-particle charge 
gap E c can be different from the optical gap E opt . 

All DMRG methods have a truncation error which is 
reduced by increasing the number m of retained density 
matrix eigenstates (for more details, see Refs. |l| and |^). 
Varying m allows one to compute physical quantities for 
different truncation errors and thus to obtain error es- 
timates on these quantities. For some quantities, espe- 
cially eigenenergies, it is possible to extrapolate the re- 
sults to the limit of vanishing truncation error and thus 
to achieve a greater accuracy. I have systematically used 
these procedures to estimate the precision of my numer- 
ical calculations and adjusted the maximal number m of 
density matrix eigenstates to reach a desired accuracy. 
This is especially important for DDMRG calculations as 
truncation errors in dynamical spectra can greatly vary 
as a function of the frequency ui for fixed m. In this work 
the largest number of density matrix eigenstates used is 
to = 600. For all numerical results presented here DMRG 
truncation errors are negligible. 

In the following three subsections I demonstrate the 
finite-size scaling technique and the accuracy of DDMRG 
on three special limits of the Pcicrls-Hubbard model. 
Then the optical conductivity of a Mott-Peierls insula- 
tor is presented and discussed in the last subsection. 



A. Peierls insulator 

For U — the Hamiltonian ( |33| ) describes a system 
of independent electrons, which can be solved exactly 
for any value of A, boundary conditions, or lattice size. 
This provides us with a perfect test case for the DDMRG 
method. I have checked that DDMRG can reproduce 
the optical spectrum of this system on lattices with sev- 
eral hundred sites, for any frequency u>, and with rela- 
tive errors as small as 10 ~ 4 using only a few hundred 
density-matrix eigenstates. This demonstrates that one 
can obtain almost exact results for the optical conductiv- 
ity of finite one-dimensional systems such as the Pcicrls- 
Hubbard model using DDMRG. 



In the thermodynamic limit the Hamiltonian (33) de- 
scribes a Peierls (band) insulating phase for A / and 
U = 0. The optical gap E opt , the charge gap E c , and the 
spin gap E s equal 2|A|. Optical excitations are made 
of one hole in the valence band and one electron in the 
conduction band. The optical conductivity is given by 




CO 



FIG. 1: Optical conductivity of a Peierls insulator with 
A = 0M and a broadening rj/t = 0.05. Both the DDMRG 
result for a 128-site chain and the exact result ( |47| ) in the 
thermodynamic limit are shown. 



for 2|A| < fiu < At and is zero elsewhere. This optical 
spectrum contains a single band of width At — 2|A| with 
square-root divergences at both band edges. These di- 
vergences are a typical feature of a one-dimensional band 
insulator. The convolution of Eq. (|47| ) with a Lorentzian 
distribution of width r]/t = 0.05 is shown in Fig. [I] for 
A = 0.6i. Both divergences are replaced by maxima at 
flu fts 2 A = 1.2t and fiu> ~ At. In Fig. [j] I also show the 
optical conductivity calculated with DDMRG on a 128- 
site lattice with the same broadening. We see that the 
finite-system optical spectrum is indistinguishable from 
the infinite-system spectrum. The broadening rj/t = 0.05 
satisfies the condition ( p6| ) and thus conceals the finite- 
size effects as discussed in Sec In this case a broad- 
ening r)(N)/t — 6.A/N is enough because the spectrum 
band width is smaller than At. 

With this size-dependent broadening rj(N) one can use 
Eq. ( p5| ) to extrapolate the finite-size DDMRG results to 
the thermodynamic limit. For instance, for fouj = 2.d>t I 
have obtained <j\(lo) = 0.245 (in units of ae /K) using 
data for systems with up to N — 256 sites [i.e., with a 
broadening down to rj(N)/t — 0.025], in excellent agree- 
ment with the exact result 0.243. If we did not know the 
exact result (^), we could nevertheless determine the 
existence of square-root divergences at both band edges 
using a scaling analysis of the maxima in the DDMRG 
spectra. For instance, the height of the low-energy maxi- 
mum (close to %uj = 1.2t) diverges as i/y/fj for r/(N) — > 
[see Fig. @(a)]. Moreover, the position of the maximum 
tends from above to the optical gap E opt = 2 A = 1.2t 
for N — > oo [see Fig. §j(b)]. As explained in Sec. IV these 



C7l (w) 



ae 2 (2A) 2 (4t) 2 



27i(M 2 \/[(M 2 " (2A) 2 ][(4<) 2 - {fiiof 



(47) 



scaling properties correspond to a square-root divergence 
at the band edge. Figure ||(b) also shows the finite- 
size optical gaps E opt (N) calculated with the DDMRG 
method for individual excited states. They tend to the 



exact result E, 



opt 



1.2t for N — > oo as expected. 
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FIG. 2: (a) Maximum <r max of the optical spectrum <Tx(u>) 
as a function of the broadening T)(N). (b) Position of the 
spectrum maximum i? ma x = TiuWx (square) and optical gap 
-Eopt (circle) as a function of system size N. In both figures 
filled symbols correspond to the Mott-Hubbard insulator [U = 
it, r/(N)N = 12. 8t] and open symbols to the Peierls insulator 
[A = 0.6t,r)(N)N = 6.4t]. 



B. Mott-Hubbard insulator 

For A = the Peierls-Hubbard model ([33|) becomes 
the one-dimensional Hubbard model at half filling. For 
U > this model describes aJMott-Hubbard insulator 
with gapless spin excitations .Ej The optical conductiv- 
ity of this system has recently bgen determined using 
DDMRG and analytical methods Here I only summa- 
rize the most important results and give more informa- 
tion about the finite-size scaling analysis carried out in 
this previous work. 

In the half-filled Hubbard model an optical excitation 
is made of a pair of spinless bosonic excitations carrying 
opposite charges in the lower (holon) and upper (dou- 
blon or anti-holon) Hubbard bands, respectively. As in a 
Peierls insulator, the optical spectrum consists of a single 
band but its width is larger, about 8t. A second distinc- 
tive feature of this spectrum is a square-root vanishing 
a\(uo) ~ y/Twj — E opt at the band threshold E opt . There 
is also a tiny peak in the middle of the band, at least for 
U > U. 

In Ref. ^ it is shown that for weak coupling (U < 3t) 
and in the strong-coupling limit (U/t — > oo) the finite- 
system optical spectra calculated with DDMRG agree 
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FIG. 3: Optical conductivity of the Hubbard model with 
U = 3t for several values of r\. Results for Tj(N) > have been 
calculated with DDMRG on iV-site lattices with jj(N)N = 
12. St. For rj = the field-theoretical result for an infinite 
system is shown. 



perfectly with the analytical results obtained in the ther- 
modynamic limit using a field-theoretical approach and 
a strong-coupling analysis, respectively. For instance, 
Fig. H shows the low-energy parts of DDMRG spectra cal- 
culated for three different lattice sizes at U = 3t and the 
corresponding field-theoretical spectrum for an infinite 
system. A size-dependent broadening ij(N)/t = 12.8/iV 
is used in this case. One clearly sees the convergence of 
the finite-size spectra toward the field-theoretical result 
as rj(N) decreases. To make a quantitative comparison 
one can calculate the convolution of the field-theoretical 
spectrum with a Lorentzian of width 7/ satisfying the con- 
dition 



IV. One finds then that 



as discussed in Sec. 
finite-size effects are completely concealed by the broad- 
ening even for relatively small system sizes. For instance, 
it is shown in Fig. 3 of Ref. || that the low-energy optical 
spectrum calculated on a 128-site lattice for U = 3t is in- 
distinguishable from the field-theoretical spectrum with 
the same broadening r\jt = 0.1. In the strong-coupling 
limit DDMRG and analytical results agree even better 
and finite-size effects are no longer visible for systems as 
small as N = 32. 

For other coupling strengths (4 < U/t < oo) it is neces- 
sary to analyze the scaling of the finite-system DDMRG 
spectra to determine the optical conductivity of the Hub- 
bard model in the thermodynamic limit. Using numerical 
results for lattices with up to N = 256 sites [i.e., a reso- 
lution rj(N)/t down to 0.05], I have found that for all U/t 
the optical conductivity at the lower band edge has the 
qualitative behavior described by Eq. (|32|): crx(u)) van- 
ishes as yjfiu) — E opt at the band threshold and there is 
a maximum in ai(uj) at a frequency u> = (1 + 7)£ , op t/?i, 
where 7 is a small number. Field theory predicts the 
same behavior with 7 1=3 0.24 in the weak-coupling limit 
while the strong-coupling analysis gives a maximum at 
hui = U = E opt + 4t, and thus 7 vanishes as t/U for 
U ^> t. The distance 7-E pt between the spectrum thresh- 
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old and the maximum increases with U/t. For U > At 
this distance is large enough to determine the finite-size 
scaling of the lower band edge using systems with up 
to N — 256 sites. As an example, Fig. ||(a) shows the 
low-energy maximum in the optical spectrum ai(u>) cal- 
culated with DDMRG for U = At as a function of r)(N). 
The contrast between the Mott-Hubbard insulator and 
the Peierls insulator is striking and the maximum in the 
Mott-Hubbard insulator optical spectrum clearly tends 
to a constant for r){N) — > 0. [For U = At the optical 
gap of the Hubbard model is comparable to that of the 
Peierls insulator with A = 0.6t, so that a direct compar- 
ison of both systems is relevant.] In Fig. ^(b) one sees 
that the finite-size optical gans calculated with DDMRG 
converge to the exact resultli-3 E op t = 1.287t in the ther- 
modynamic limit, but the maximum tends to a higher 
energy Tllo w 1.7t. Therefore, one can conclude that 
there is no divergence at the optical conductivity thresh- 
old hto — E opt . Moreover, it is possible to confirm that 
<7i(tL>) vanishes as y/hui — E opt at the lower band edge us- 
ing either a similar scaling analysis for the derivative of 
DDMRG spectra or a direct comparison with the convo- 
lution of functions-such as Eq. (|32| ) or the field-theoretical 
optical spectrum. II 

For very weak coupling one would need to calculate 
cti(w) for very large system sizes in order to perform 
the same scaling analysis. Because .Eopt vanishes ex- 
ponentially with U/t and the scaling analysis must be 
performed in the asymptotic regime rj(N) < jE opt , the 
required system sizes N increases exponentially as t/E op t 
for U/t — ► 0. Fortunately, it is not necessary to carry out 
this analysis for the Hubbard model because the optical 
conductivity of the weak-coupling field theory is already 
in excellent agreement with the optical conductivity of 
the lattice model for U = 3t. 
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FIG. 4: Reduced optical conductivity C,[cJa"i(u;)] as a func- 
tion of hfl = huj — U in the strong-coupling limit U ^> t for 
A = 0M. The solid line is the DDMRG result, Eq. @, for a 
128-site lattice with a broadening r\/t = 0.1. The dot-dashed 
line is a Lorentzian distribution of width r\/t = 0.1 centered 
at hw = U. The two dashed lines represent the analytical 
result @ for the continuum of an infinite system (77 = 0). 
Note the logarithmic scale of the ordinate axis. 



decreases exponentially with increasing distance. Thus, 
this strong-coupling limit of the Peierls-Hubbard model 
is different from the two limiting cases discussed previ- 
ously and from the general case presented in the next 
section. 

In the thermodynamic limit the optical conductivity 
can be calculated analytically using some reasonable as- 
sumptions .E_3 If < |A| < 2t, the spectrum consists of 
two bands for 2|A| < \hcj -U\<2t 



g e 2 a y/[(nny - (2Ay][(At) 2 - (nnf 



hu>\h£i\ 



(48) 



C. Strong-coupling limit 

In this section I discuss the special case of a Mott- 
Peierls insulator (A 7^ 0, U > 0) in the strong-coupling 
limit U/t — > 00, for which the. shape of the optical spec- 
trum is known analyticallyJlj In this limit there is ex- 
actly one electron on each site in the ground state of the 
Peierls-Hubbard Hamiltonian ([33]) . An optical excitation 
moves one electron from a site to another and thus creates 
a double occupation (doublon) and an empty site (holon). 
Therefore, the optical gap is of the order of U. These el- 
ementary charge excitations are spinless bosons as in the 
Hubbard model. The properties of the spin degrees of 
freedom are determined by an effective Heisenberg model 
with alternating exchange couplings J\ ~ (t + A/2) 2 /U 
and J2 ~ (t — A/2) 2 /U. The spin gap E s vanishes in 
the limit U — > 00. However, as there is a gap in the spin 
excitation spectrum for any finite U ft (see also next sec- 
tion), the structure of the spin ground state in the limit 
U/t — > 00 is actually similar to that of a gapped state. 
For instance, the antiferromagnetic spin-spin correlations 



where hfl — Hlu — U, and a ^-peak at Tiu = U 

7r 5Tr e 2 ai 2 
0-1(0; ) = — — — o(huj - U) 
nU 



(49) 



in the middle of the gap 4|A| separating the bands. For 
A — > one recovers the optical spectrum of the Hubbard 
model in the strong coupling limit, which consists ofr-a. 
single band and a <5-peak in the middle of this band]§t3 
The prefactors go and are spin form factors given by 
ground-state spin correlation functions. They are func- 
tions of the effective exchange coupling ratio J2/J1 and 
thus of (5 = |A/2i|. Assuming a dimer spin ground state 
(i.e., J\ > and J2 = 0) one obtains go = 9/4 and 



1 + 32(5 + 62<5 2 + 32<5 3 + S 4 
4(1 + 5) 2 



(50) 



This result becomes exact in the dimer limit |A| = 2t, 
where g n — 8. For A — ► the dimer spin ground state 
does not give the correct form factors because it is known 
exactly that go + g v = A ln(2) and it was found numeri- 
cally that g^/go ~ 10~ 2 (see Ref. H). 
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Figure |] shows the reduced optical conductivity 
C v [u)ax{ijj)\ calculated using DDMRG on a 128-site lat- 
tice for A = 0.6t and 37/t = 0.1. A logarithmic scale 
is used to make visible the weak bands on both sides of 
the strong central peak. In Fig. ^ one can recognize the 
spectral shape predicted by the strong-coupling analysis. 
To make a quantitative comparison, however, it is first 
necessary to determine go and g n using the finite-size 



scaling analysis of Sec. IV. Here I use a size-dependent 
broadening rj(N)/t = 12.8/iV as for the Hubbard model 
because the spectral width is also of the order of 8t. For 
hcj = U - 2t DDMRG results for C v [ojui(uj)] tend to 
0.78 (in units of e 2 at/h 2 ) for N — ► 00. Comparison with 
Eq. ( ^8| ) then yields 50 ~ 2.2. In Fig. || I also show 
the two bands (^) with this value of go (without broad- 
ening). The agreement with the finite-system DDMRG 
spectrum is excellent. The small deviations visible close 
to the band edges are due to the different broadening used 
for the numerical result (r\jt = 0.1) and for the analyti- 
cal result (77 = 0). They vanish if the same broadening is 
used in both calculations. Once more this confirms that 



a broadening satisfying Eq. (26) hides most finite-size ef- 
fects in this model as already shown by other examples 
in Sec. VA of this paper and in Refs. || and ^. In the 



DDMRG spectra C v [u>cri(u>)] the height of the central 
peak diverges as 2A2tjr\ (in units of e 2 at/Ti 2 ) for increas- 
ing rj(N) but its position does not change. This confirms 
that it corresponds to a J-peak at huj — U and gives 
an estimate </v ~ 2.42. This 5-peak broadened with a 
Lorentzian distribution of width rj/t — 0.1 is also shown 
in Fig. |[ One sees that the agreement with the DDMRG 
result is perfect. A similar finite-size scaling was per- 
formed-.to determine the form factor g w in the Hubbard 
model.B 

One notes the surprisingly good agreement between 
the form factors determined numerically with DDMRG 
(go ps 2.2 and g v ss 2.42) and those obtained using the 
approximation of a dimer spin ground state (go = 2.25 
and g^ w 2.52). For the value A = 0.6t used in this 
example, the ratio J2/J1 ~ 0.29 of the effective exchange 
coupling is already quite small and thus the dimer spin 
ground state is probably a very good approximation of 
the actual spin ground state. 



D. Mott-Peierls insulator 

The optical conductivity of the Peierls-Hubbard model 
is not known for general interaction parameters U > 
and 2t > |A| > 0. In thia.regime the system is in a Mott- 
Peierls insulating phaseEj: both a periodic lattice poten- 
tial (i.e., the alternating hopping terms) and electronic 
correlations contribute to the formation of a charge gap 
E c > and there is a finite spin gap E s > 0. Numerical 
investigations of the charge and spin gaps and of static 
correlation functions reveal no phase transition at finite 
U and intermediate A (see also Ref. |l(^ and references 
therein). Thus the entire parameter space (0 < U/t < 00, 




FIG. 5: Charge (diamond) and spin (circle) gaps extrapo- 
lated to the thermodynamic limit, (a) As a function of U for 
A = OM. (b) As a function of A for 17 = 41 



< |A| < 2i) belongs to a single Mott-Peierls insulating 
phase. Figure |^ shows charge and spin gaps as a function 
of U and A. These gaps have been calculated on lattices 
with up to N = 400 sites using DMRG and extrapo- 
lated to the thermodynamic limit. The charge gap of the 
Mott-Peierls insulator is always larger than the gap of 
the Mott-Hubbard and Peierls insulators in the A = 
and U = limits, respectively. The spin gap is always 
smaller than the charge gap in the Mott-Peierls phase. 

In the thermodynamic limit the optical gap E opt is 
equal to the charge gap. The nature of the optical ex- 
citations in the Mott-Peierls insulator is not well under- 
stood. Despite the obvious difference between charge and 
spin excitation energies, E s < E c , it is not even known 
if there is a spin-charge separation for single-particle ex- 
citations. Optical excitations could consist of a pair of 
fermionic quasi-particles with opposite spins ±<r a nd op - 
posite charges ±e as in a Peierls insulator (Sec. |V A| ). 
They could as well be made of two spinless bosonic ex- 
citations carrying oppos ite c harges ±e as in the Mott- 
Hubbard i nsula tor (Sec. VB) and in the strong-coupling 
limit (Sec. |vtj ). 

The investigation of spin and charge gaps and static 
correlation functions clearly shows that the three spe- 
cial cases described in the previous sections are singular 
limits of the Peierls-Hubbard model. Unsurprisingly, I 
have found that the optical conductivity in the Mott- 
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FIG. 6: Reduced optical conductivity C v [u)ai(u>)] calculated 
with DDMRG [see Eq. @] on a 128-site lattice (-q/t = 0.1) 
in the strong-coupling regime for A = 0.6t as a function of 
HQ = Tilu — U. Note the logarithmic scale of the ordinate axis. 



2.0 




FIG. 7: Optical conductivity a\{uj) calculated with DDMRG 
on a 64-site lattice (rj/t = 0.2) for A = 0.6t and several values 
of U. 



Peierls phase is unlike the simple spectrum observed in 
these limits. [All optical spectra presented in this section 
have been calculated using DDMRG and the finite-size 
scaling analysis has always been performed using a size- 
dependent broadening rj(N)/t = 12.8/JV.] 

For large but finite U the optical spectrum consists 
of three bands: a narrow band with a strong singularity 
around hu) — U and one weak band on each side of this 
central peak. The singularity seems to be made of two 
very close power-law divergences which merge to form 
the single isolated <5-peak ( J49| ) in the U/t — > oo limit. 
The optical spectrum starts as y^hu — -E"opt a t the lower 
band edge E opt for all |A| < 2t. Figure ^ shows the re- 
duced optical conductivity C 7? [o;iti(w)] for U — AOt and 
A = 0.6t. The spectrum looks very similar to the spec- 
trum for U/t — > oo, which is again shown in this figure. A 
finite-size-scaling analysis shows however that the strong 
central peak is not a (5-function for U = AOt but a nar- 
row band with at least one singularity diverging as 7/~ ' 8 . 
The spectra in Fig. ^ are made of three bands: the gaps 
between the bands appear as local minima on both sides 
of the central peak because of the relatively large broad- 
ening used {rj/t — 0.1). The finite-size scaling analysis 
confirms the existence of these gaps. For decreasing pa- 
rameters U or |A| first the lower gap, then the upper 
gap close. Therefore, the number of bands in the opti- 
cal spectrum of the Mott-Peierls insulator is not constant 
but depend on the interaction parameters U and A. 

The evolution of the optical conductivity as a function 
of U is very interesting. For decreasing U ft one observes 
that the central peak breaks into two peaks appearing 
as local maxima in the broadened spectrum of finite- 
size systems. The first peak (at the lowest energy) takes 
over most of the optical weight of the central peak. Its 
weight decreases progressively with decreasing U/t but 
remains strong even for small U. In Fig. |7, it is clearly 
visible (at Tiuj > At) even for U = 2t (with A = 0.6i). 
This peak corresponds to a power-law divergence within 



a band with an exponent that tends to —1/2 for U — ► 0. 
The peak position moves to lower energy as U decreases 
and reaches hcu — At for U = 0. Therefore, the central 
peak observed at strong coupling U ^> t corresponds to 
the upper square-root divergence in the Peierls insulator 
spectrum ([i"7|). (In Fig. [7] this divergence is barely visible 
as a local maximum at hut ss At because of the relatively 
large broadening rj/t = 0.2 used.) The second peak has 
very little optical weight (it is not visible in Fig. [7]) and 
I have not been able to determine its structure. 

For U 3> t the optical weight is distributed symmet- 
rically around the central peak, as seen in Fig. ||. As U 
decreases, there is a progressive transfer of optical weight 
from high frequency (above the central peak) to low fre- 
quency (below the central peak), see Fig. 0. The high- 
frequency spectrum becomes very weak for small U ft but 
it completely disappears only at U = 0. It is difficult 
to determine the spectral width for general parameters 
because the optical conductivity is very weak and van- 
ishes smoothly at high frequency. I estimate that the 
width of the spectrum lies between At and 8t for U > 0. 
The smallest width is reached for large A and small U, 
the largest for small A and large U. The low- frequency 
spectrum becomes stronger as U diminishes. The local 
maximum below the central peak (seen in Fig. ^) pro- 
gressively rises, moves closer to the lower band edge, and 
transforms into a strong narrow peak, visible in the spec- 
tra shown in Fig. (at ficu < At). For small enough U/t 
this low-energy peak contains more optical weight than 
the central peak. For U — > the low energy peak be- 
comes the square-root divergence of the Peierls insulator 
spectrum (47) at the band threshold E opt . 

For U > 0, however, my results suggest that the op- 
tical spectrum always vanishes smoothly at the optical 
gap. I think that the low-energy optical spectrum of 
the Peierls-Hubbard model at weak coupling has a qual- 
itative behavior similar to that of the Hubbard model: 
as hio — E opt decreases, <j\ (u>) first appears to diverge 
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FIG. 8: Optical conductivity ujC v [ct((.u)] calculated with 
DDMRG [see Eq. @] on a 128-site lattice {rj/t = 0.1) in 
the small gap regime: Mott-Hubbard insulator with U = 3t 
(dashed), Peierls insulator with A = 0.3i (dot-dashed), and 
Peierls-Hubbard insulator with U = 2.3t and A = 0.15i 
(solid). The optical gaps are E op t = 0.631£, 0.6i, and 0.704t, 
respectively. 



as (hu> — -E opt ) -1 / 2 , then goes through a maximum just 
above the optical gap E opt , and vanishes smoothly for 
Huj — ► E opi . For large enough U it is possible to carry 
out a finite-size scaling analysis similar to th e one per- 
forms for the Hubbard model (see Sec. VB ). Thus it 
is possible to check that the low-energy spectrum maxi- 
mum is finite and lies at a higher energy than the optical 
gap, and to show explicitly that <j\{ui) ~ yjfioj 



Huj — E. 



opt 



0+. 



-Eopt for 



For smaller U it becomes increasingly difficult to dis- 
tinguish a smooth spectrum with a truncated divergence 
from a true divergence. For instance, Fig. @ shows 
the low- frequency optical conductivity wC,[(Ti(u)/w] for 
U = 2.3i and A = 0.15i with a broadening rj/t = 0.1 
(N = 128 sites). For comparison, I also show the spec- 
tra in two limits discussed previously, the Mott-Hubbard 
insulator [o\{uj) 
tor [(7i(u;) 



yjhuj — Eopt ] an d the Peierls insula- 
with similar optical gaps 



UJ 



E, 



opt 



{E pt/t = 0.6 — 0.7) and the same broadening 77. Clearly, 
the Mott-Peierls insulator spectrum looks like an inter- 
mediate case between the spectra observed in both lim- 
iting cases. The position of the maximum in the Mott- 
Peierls insulator spectrum tends to 0.78i for rj(N) — > 
while the charge gap (and optical gap) equals 0.704?; in 
the thermodynamic limit. Certainly, there is no diver- 
gence in the low-energy optical spectrum. However, the 
maximum seems to diverge as 1/ ^fr\ even for the small- 
est broadening I have used {rj/t — 0.05). Thus this spec- 
trum seems to be qualitatively similar to a function such 
as Eq. (|32j), but the maximum is so close to the optical 
gap that broadenings 77 significantly smaller than 0.05/; 
(i.e., system sizes much larger than N = 256) would be 
necessary to reach the asymptotic regime as discussed in 
Sec. [V. For the same reason it is not possible to deter- 



FIG. 9: Optical conductivity eri(w) of a Mott-Peierls insula- 
tor calculated with DDMRG on a 128-site lattice {rj/t = 0.1) 
for U = 2.3t and A = 0.15t (dashed) and the fitted field- 
theoretical spectrum for Mott insulators with f3 2 = 0.58 and 
the same broadening rj/t (solid). 



mine how the spectrum vanishes for hu> — £opt — * + in 
such a case 

In the Hubbard model it is possible to confirm the 
absence of a singularity and the square-root vanishing at 
the band threshold even if the optical gap is as small as 
-Eopt = OAt, because we know the optical spectrum of 
an infinite system for E opt — > from field theory.Q The 
field theory approach does not only apply to the Hubbard 
model, but more generally, gives the low-energy optical 
spectrum o£pae-dimensional Mott insulators with small 
Mott gaps.BEj The different spectral functions depend 
only on an interaction parameter /3 2 < 1. In addition, 
the optical gap E opt > and a normalization constant 
set the frequency scale and the conductivity scale. For 
1/2 < /3 2 < 1 these optical spectra described truncated 
square-root divergence with a square-root vanishing at 
the band threshold as in Eq. @. For 1 = 1 (Hubbard 
model) the spectrum has the shape shown in Fig. || with a 
maximum at 1.24-B opt . As 1 decreases the peak becomes 
sharper and the maximum moves closer to the band edge. 
For P 2 = 1/2 the optical spectrum is similar to that of 
a Peierls band insulator with a square-root divergence at 
the band threshold. 

Therefore, this field theoryH'EJ can describe the op- 
tical spectrum in both limiting cases (Mott-Hubbard 
and Peierls insulators) of the Peierls-Hubbard model in 
the small gap regime, and the field-theoretical spectrum 
evolves continuously from one limit to the other with f3 2 
going from 1 to 1/2. Using 1 > 1 > 1/2, the optical 
gap -E opt , and the normalization constant as fit parame- 
ters, I have compared the low-energy optical conductivity 
calculated with DDMRG for small gaps (E opt < 0.71t) 
to field-theoretical spectra with similar broadening as 
explained in Sec. [lV|. For instance, I show in Fig. ^ 
the DDMRG spectrum for the lattice model (|3^) with 
U = 2.3t and A = 0.15t and the fitted field-theoretical 



15 



spectrum with 1 
huj = 1.2* w 1.7E, 



0.58. Both spectra agrees up to 
opt . Generally, I have found that the 
optical spectrum of a Mott-Peierls insulator can be fit- 
ted by a field-theoretical spectrum with /3 2 > 1/2 over 
a finite frequency range, from oj = to a frequency w 
which lies between the position of the low-energy maxi- 
mum and 2E opt /h. (Naturally, for U — > the best fit is 
always obtained with (3 2 — 1/2.) Therefore, I think that 
for any U > (and < |A| < 2t) the optical spectrum 

0+. 



E opt for hu> 



E, 



opt 



vanishes as yffou 

Note that I do not assume that the field-theoretical cal- 
culations in Refs. ^| and 16 are also valid for the Peierls- 
Hubbard model with general interaction parameters. Ac- 
tually, there are visible discrepancies starting at rather 
low energy between field theory and DDMRG results for 
the lattice model, as shown in Fig. ^. The agreement be- 
tween field theory and DDMRG results in the region of 
the band threshold simply means that the optical spectra 
in Mott-Pcicrls insulators and in one-dimensional Mott 
insulators have similar shapes just above the optical gap. 

Finally, it is interesting to examine the evolution of 
the optical spectrum from weak to strong bond alterna- 
tion for fixed U. It has been shown in Ref. || (see also 
Sec. VB) that for A = (Hubbard model) the spectrum 
consists of a single band with a maximum close to the 
lower band edge and a tiny peak in the center (at least 
for U > At). If |A| increases one observes in Fig. |Io|(a) 
that the maximum moves closer to the optical gap and 
corresponds to a sharper peak. The optical spectrum 
still starts as <J%oj — E opt at the band threshold as dis- 
cussed above. The central peak, which is too weak to be 
seen in Fig. |l0|(a) for A = 0, becomes rapidly stronger 
as |A| increases and is clearly visible for A = Q.At. As 
discussed previously this peak becomes a 5-function in 
the strong-coupling limit U t and corresponds to the 
upper square-root divergence of the Peierls insulator if U 
vanishes. For moderate |A| the ratio between the hop- 
ping integrals r(A) = (t - |A|/2)/(f + |A|/2) is not too 
small and the optical weight is mostly concentrated below 
the central peak. If this ratio becomes small, however, 
the central peak becomes the spectrum dominant fea- 
ture, see Fig. 10Kb) . The proportion of the optical weight 



which is in the central peak increases as 1 — r (A) for 
r(A) — > 0. A finite-size scaling analysis confirms however 
that this peak is not a (5-function but is still a power-law 
divergence within an excitation band. Only in the dimer 
limit |A| = 2t [r(A) = 0] the optical spectrum is made 
of a single <5-peak, which corresponds to the excitation of 
Frenkel excitons localized on a dimer. 

In summary, the optical spectrum of a Mott-Peierls 
insulator consists of one or more bands with a total spec- 
tral width ranging from At to 8t. The distinctive features 
of the spectrum are a square-root vanishing of <j\ (uj) at 
the lower band edge and a peak due to a power-law sin- 
gularity around the middle of the spectrum. For strong 
couplings [U 3> t and r(|A|) <C 1] most of the optical 
weight is in the central peak, while for weak couplings 
[U <C At and r(| A|) > 1/2] it is mostly concentrated in a 
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FIG. 10: Optical conductivity ai(iv) calculated with 
DDMRG on a 64-site lattice (rj/t = 0.2) for U = At and 
several values of A. (a) For r(A) > 0.5. (b) For r(A) < 0.5. 



narrow peak just above the optical gap. In the limit of 
a vanishing gap (U — > and A — > 0) this narrow peak 
becomes a Drude peak at u = 0. For intermediate cou- 
plings most of the optical weight is distributed over a 
broad frequency range between the optical gap and the 
central peak and <j\{lo) goes through a maximum in this 
range. 

The central peak always appears at an energy larger 
than the bare band width At. For_Darameters which are 
realistic for conjugated polymers,EB most of the optical 
weight lies below this peak. Therefore, I think that it 
is not possible to observe such a structure in the optical 
spectrum of conjugated polymers, because it occurs at a 
too high energy (> 10 eV) and its intensity is too weak. 

The optical spectrum of Mott-Peierls insulators is un- 
like that of Mott-Hubbard and Peierls insulators. How- 
ever, most of its main features are found in the strong- 
coupling limit investigated in Ref. [TBI and discussed in 



Sec. VC. This suggests that the optical excitations of 
a Mott-Peierls insulator could be made of a pair a spin- 
less bosonic excitations with opposite charges as in the 
strong-coupling limit (and in a Mott-Hubbard insulator) . 
Nevertheless, understanding the nature of the system el- 
ementary excitations requires the study of additional dy- 
namical properties, such as the one-particle Green's func- 
tions. The DDMRG method will enable us to carry out 
this further investigation. 
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VI. CONCLUSION 

In this paper I have presented a dynamical DMRG 
method which allows one to calculate the optical con- 
ductivity of one-dimensional correlated electron systems 
on large lattices with great accuracy. The DDMRG ap- 
proach to the calculation of dynamical properties is es- 
sentially an application of the standard DMRG algorithm 
for ground-state calculations. Therefore, both methods 
have the same advantages but also the same limitations. 
In particular, DDMRG will directly benefit from recent 
and future improvements of DMRG such as the use of 
additional symmetries. 

With DDMRG it is possible to calculate the dynami- 
cal response of correlated systems with hundreds of sites 
and particles. Relative errors of the order of 10 -4 can 
be achieved for the optical spectrum of finite systems 
with N ~ 10 2 sites. Using a finite-size scaling analy- 
sis based on a size-dependent broadening of the discrete 
finite-system spectra, one can then calculate a dynamical 
spectrum in the thermodynamic limit with a resolution 
of the order of 1% of the spectral width and even inves- 
tigate singularities in a continuum. 

The DDMRG approach can be used for various dynam- 
ical quantities, such as dynamical spin-spin correlation 
functions or single-particle Green's functions. It can also 
be applied to other lattice quantum many-body models, 
in higher dimension, including spin or boson degrees of 
freedom, or long-range interactions. The correction vec- 
tor DMRG method has been used to calculate non-linear 
dynamic response functions, such as third-order dynam- 
ical polarizabilitiesE Similarly, the variational principle 
presented in Sec. II can be generalized to dynamical cor- 
relation functions describing these non-linear responses. 
Thus I believe that it is possible to develop an efficient 
DDMRG method for calculating these quantities. A lim- 
itation of the DDMRG approach is the restriction to zero 
temperature. It would be desirable to extend the vari- 



ational principle and the DDMRG approach to finite- 
temperature dynamical properties. 

The computational resources used by the DDMRG 
method are relatively modest. For instance, all calcula- 
tions presented in this paper were carried out on worksta- 
tions with a single 500 MHz Alpha processor and required 
less than 1 GByte of memory. It would be easy and very 
efficient to run DDMRG on a parallel computer as calcu- 
lations for different frequencies are almost independent. 
This would permit one to investigate much larger or more 
complicated systems than in this work. 

In summary, the DDMRG method and the finite-size 
scaling technique for dynamical spectra appear extremely 
accurate and versatile. They provide a powerful ap- 
proach for investigating the dynamical properties in low- 
dimensional lattice quantum many-body systems. 

Finally, it should be kept in mind that the variational 
principle for dynamical correlation functions and their re- 
lated excited states is completely independent from the 
DMRG method. Therefore, it is possible to combine this 
principle with other variational methods to calculate dy- 
namical properties. For instance, one could build a trial 
wavefunction \ip({Xi})) — i£({Aj})|^>o)i where R({Xi}) is 
an operator depending on a few parameters Xi , such that 
the calculation of lF"({Aj}) = Wyt j7) (w, ^({-^i})) reduces 
to the evaluation of ground-state correlation functions. 
Then the minimization of W^({Aj}) with respect to the 
variational parameters {A^} would give a lower bound 
and an approximate value of the dynamical spectrum 
I a (uj + iif). 
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